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Abstract 



A method for the evaluation of the e-expansion of multi-loop massless Feynman in- 
CN ' tegrals is introduced. This method is based on the Gegenbauer polynomial technique 

and the expansion of the Gamma function in terms of harmonic sums. Algorithms 
f-^. , for the evaluation of nested and harmonic sums are used to reduce the expressions 

to get analytical or numerical results for the expansion coefficients. Methods to 
increase the precision of numerical results are discussed. 
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1 Introduction 



One of the most important aims of particle physics is the test of models, 
particularly the standard model. Testing the standard model means, on one 
hand, to increase the accuracy of its parameters and, on the other hand, the 
search for deviations of its predictions from measured quantities to discover 
new effects. To do this experiments are accomplished with increasing accuracy 
of measurement. Therefore theoretical calculations have to be performed with 
increasing numerical precision. 

A parameter of the standard model of particular interest is the strong coupling 
constant a s . It can, for example, be determined via the so-called i?-ratio, the 
normalized total cross section of the process e + e~ — > hadrons, 

. . o~(e + e~ — > hadrons) . , 

R(s) = -t+ 1 

a(e + e — > //+ /x ) 
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In order to achieve high precision one has, within perturbative quantum field 
theory, to calculate high orders of the perturbation series. This means, we need 
to calculate large numbers of more and more complicated Feynman diagrams. 
Therefore different methods for the evaluation of such multiloop Feynman 
diagrams have been developed. (See e.g. [1,2,3].) 

The quantity R(s) is nowadays known to 0{ot z s ) precision [4,5] and first results 
in O(aj) have been published [6]. For a review see Refs. [7,8]. 

A prescription to solve problems with a large number of complicated Feynman 
diagrams is to reduce them to a set of master integrals, e.g. by means of the 
traditional integration by parts method, and then evaluate these. 

The aim of the project described in this article is to develop a method for 
the calculation of massless integrals. More precisely, we want to consider the 
planar (P) and non-planar (N) three-loop propagator diagrams: 




They can be considered as the master integrals of the corresponding class of 
Feynman diagrams, which are "building blocks" of the calculations for R(s) 
mO(aj). 

We are interested in the ^-expansion of the diagrams, where e is the parameter 
of dimensional regularization. This means, the integrals are formally calculated 
in D = 4 — 2e dimensions and then expanded in e. By means of this procedure 
UV divergencies of diagrams manifest themselves as poles in e and can then 
be cancelled by means of renormalization. 

It is thus essential for renormalization to know the pole structure of a Feynman 
diagram and the finite part of the expansion. However, when one has to deal 
with products of graphs with poles as in the case of the following diagram, 




which appears within the mentioned calculations of R(s), one has to calculate 
higher e orders of the second factor in order to ensure, that the result is correct 
up to the finite term. This is why we are interested in higher orders of the e 
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expansion of the diagrams under consideration. 

The 0(e v ) terms of the diagrams iV and P are necessary for the computation 
of the (3 function of the scalar 4 model at five-loop level. They have been 
known since long from Ref. [9] and [10]. 

The calculation of the planar diagram is somewhat simpler than that of the 
nonplanar one, because it could be related to the generic massless two-loop 
two-point function by means of integration by parts [11]. The latter was stud- 
ied in many works, e.g. [10,12,13,14,15,16] and references therein. Recently, in 
Ref. [17] it was in a sense computed analytically. Thus one can say that the 
planar diagram is known to arbitrary order and certainly to sufficient precision 
for practical needs. 

The known results for the diagrams P and N read [10,18] 

N = j^- 6 G 3 (e)(k 2 r 2 ^N s , (2) 

N s = 20 C 5 + e (68 C 3 2 - 80 ( 5 + 50 Ce) + 0{e 2 ) (3) 

and 

P = ^- ? G\e){kr 2 ^Ps, (4) 

P 5 = 20C5 + £(44C 3 2 - 80C5 + 50C 6 ) 

+ e 2 (-176CI + 132C 3 C 4 + 80C 5 - 200C 6 + 317C 7 ) + 0(e 3 ). (5) 



In the following we denote the expansion coefficients of order n by \ e n, e.g. 
N s \ e i =68C 3 2 - 80C 5 + 50Ce. 

These numbers have been used to check the method and to examine the prop- 
erties of the numerical algorithms described later. 

The determination of the e 2 coefficient of N is one of the aims of the project 
presented in this paper. An independent evaluation can be found in Ref. [19], 
the result was given in Ref. [20]. A numerical approach to the calculation of 
integrals of this type is given e.g. in Ref. [21]. 

The paper is organized as follows: Section 2 shows how a representation of the 
diagrams N and P in terms of sums can be found by means of the Gegenbauer 
polynomial method. Section 3 gives an overview over the Harmonic Sums and 
algorithms to deal with them. It is shown how the summation algorithms are 
applied to evaluate the sums representing the diagrams. For the N diagram an 
analytical solution by means of the described algorithms is not yet possible. 
In this case one can get numerical results by summing up the series. Section 4 
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reviews methods to accelerate the convergence of such processes to get better 
results. By means of these numerical methods it is possible to calculate the 
expansion coefficients with high precision. For the non-planar diagram this 
is described in Section 5. The planar topology can be calculated analytically 
with the help of the summation algorithms, which is demonstrated in Section 
6. 



2 The Gegenbauer Polynomial x-Space Technique 



The Gegenbauer Polynomial x-Space Technique (GPXT) to evaluate Feynman 
Diagrams in coordinate space was introduced in Ref. [12]. Applications can 
be found e.g. in Refs. [22,23], an extension of the technique in Ref. [24]. In 
this Section we will apply the method to the diagram N. The results for the 
planar diagram can be obtained analogously. Labeling the momenta as shown 
in Fig. 1 the integral representation of iV reads 



5 k + q + p - r 4 




1 r 2 

Fig. 1. Labeling of the momenta and vertices for diagram N 



N = wr I dD v dD « dDr 



' ■ (6) 



p2 r 2 g2 (j^ _|_ qY (j^ _j_ g _|_ p _ r y ^ _|_ py ( r _ qy ( r _ p^ 



By means of Fourier transformation 

J_ L(A) r e 2ik *d D x 

¥ ~ J {x 2 ) x ' ( ' 

where A = 1 — s = ^ — 1, the integrand is transformed to position space. The 
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integrals over the momenta can then be evaluated to delta functions: 

J e 2tp *d D p = 7r D 5(x). 



(8) 



After the evaluation of the delta functions and choosing an appropriate coor- 
dinate system (see Fig. 1, where the root vertex is denoted by 0) the denomi- 
nators get a simple and uniform shape: 



N = 26 ^ ( ^ (A+1) / d D X! d D x 2 d D x 3 d D Xi d D x 5 e~ 2 ^^ 



™2A 
x 5 



X 



(x 2 - Xi) 2X (x 3 - x 2 ) 2X (x 4 - x 3 ) 2X (x 5 - x 4 ) 2X (xi - x 4 ) 2X (x 5 - X 2 ) 2X 

(9) 



Spherical coordinates are introduced by means of 



1 2n x+1 

dDx =2 T(\ + l) rXdrdk > ( 10 ) 



where 

r = x 2 and x = —=. 



Now the propagator terms and the exponential function can be expanded in 
Gegenbauer polynomials using formulae (A. 8) and (A. 11) of appendix A. This 
leads to a representation of N as a multiple sum 

*= a «££L eeeeeee (») 

/y V / j I m n s t u 

x J dr±dxi J dr 2 dx. 2 J dr 3 dit 3 J dr 4 dx 4 J dr 5 dx. 5 
x r(A) V (j + A) C x (k • x 3 ) (k 2 r 3 y/ 2 j J+x (k 2 r 3 ) 

7" ^ 7" g 



x C x (St 3 • x 4 ) C A (x 4 • x 5 ) C A ( Xl • x 4 ) C x m (x 2 • x 5 ) 

x M{(r l5 r 2 ) M A "(r 2 , r 3 ) M A "(r 3 , r 4 ) M A m (r 4 , r 5 ) M A (n, r 4 ) M A m (r 2 , r 5 ), 



where 



v n/2 / \ n/2 

n \ I ■ ( r i r 3 



M A "(r 4 , r,) = ? - = , min A - • (12) 

max(rj, rj) A \rj / max(rj, rj) A \rj r { J 
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The angular integrations can be evaluated using the properties of the Gegen- 
bauer polynomials. Each application of the orthogonality relation (A. 4) re- 
duces two C x to one and produces a Kronecker delta. Finally one arrives at 
the expression 



"*4 / \ \ oo oo l+m 



« = J^tt E (13) 

v ' i=0 rn=0 n=\l— m\ 

l+m+n=2g , gd 



T(n + 2A) D\(l,m;n) . . 2 r\r\r\r\r\ 

X n!T(2A) (Z + A) (m + A) (n + A) 2 Ja( ^ ^ 
x M A (n, r 2 ) M A "(r 2 , r 3 ) M A > 3 , r 4 ) M A m (r 4 , r 5 ) M A (n, r 4 ) M A m (r 2 , r 5 ). 



The complicated structure of the triple sum is introduced by Eq. (A. 14), which 
is applied when there are two polynomials with the same arguments. 

To handle the radial integrals the region of integration is distributed in 120 
regions according to the relative size of ri,r 2 ,r 3 ,r^,r 5 . In each region the ex- 
pressions M A , Eq. (12), evaluate to rational expressions. For simplification k 2 
is set to 1; the /c-dependence can later be restored by dimensional arguments. 
Thus we are lead to integrals like 



rrz m+n—l rV2 rco roo poo m—n—l—2\ 

/ dr 2 r 2 2 / dr l r[ \ dr 3 j x (r 3 ) / dr 5 r^ m ^ 2X / rfr 4 r 4 2 

JO JO JO Jr 4 Jr 3 

4 r(5 - 3A) 

" (I + 1) (1 + m + n + 4) (m + 2A - 1) (1 + m + n + 6A - 4) T(4A - 4) ' 



Actually, due to the symmetry of the diagram, only 30 of these 120 terms are 
different. The calculation of these can be automated. * 

After performing all integrations and rewriting A = 1 — e the integral N is 
expressed in terms of a triple sum: 

N = C ■ N s , (14) 



We thank K. Chetyrkin for his Mathematica implementation of this procedure. 
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oo oo l+m /i \ 

iV S = EE E R(e,l,m,n) 



x 



where 



n™ n , i ' ' ' r(2-2e)r 2 (l -e) 

m=0 n=|Z-m| v / V / 

1 T(g + 2-2e) 

(I + 1 - e) (m + 1 - e) (n + 1 - e) r(# + 2 - e) 

£(g - I + 1 - g) T(g - m + 1 - e) T(g - n + 1 - e) 

r(g-l + l)r(g-m + l)r(g-n + l) ' 1 °> 



c _ r 4 (l-e)r(2 + 3e) 



R(e, I, m, n) = 



(4 7r )3(2- e )( 1 _ e )4 r (_ 4e )' 



4 



(/ + 1) (/ + m + n + 4) (m - 2e + 1) (/ + m + n - 6e + 2) 

8 

+ (I + 1) (/ + m + n + 4) (/ + m + n - Ae + 2) (/ + m + n - 6e + 2) 

+ . . . (118 similar terms). 



Finally, according to Ref. [12], this result can be cast into the so-called G- 
Form. Rewriting 

C = j±y 6 G%e) ■ (-As - 4, 2 + 32e 3 + 0(e*)) (16) 



one obtains 



N = j^G*(e)(kY 2 - 3£ N s , (17) 



N s = (_4 £ _4 e 2 + 32e 3 + C(e 4 ))A^ 5 . (18) 

N,s can now be expanded in powers of e. For each order in e a triple sum has 
to be evaluated. The calculation of these sums is the aim of the rest of the 
paper. 



3 Summation Algorithms 



In this Section the Harmonic Sums are introduced and their relationship to 
the problem under consideration is shown. 
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3.1 Harmonic Sums 



Harmonic Sums are defined as [25,26,27] 



S ™(n) = it^ (19) 



Am 
i=l 6 



which, for n = oo, leads to Riemann's zeta function: 

5 m (oo) = ((m) (m > 2). (20) 
Higher Harmonic Sums are defined [25] recursively via the relation 

n n 

•S„: ,» £^ (21) 

i=l ^ 

The number of indices of a harmonic sum is called its depth, the sum of the 
absolute values of the indices its weight. 

A negative index denotes an alternating sum, t 



n C-lV 

^) = E^ (22) 
i=i 4 

In Ref. [25] different algorithms are described to deal with expressions of the 
form 

SaM(n) = t^ a ±j b t^ (23) 
i=i ? j=i J fc=i K 

and sums over such. These are implemented in the program package SUMMER 

written in FORM [28] which is available under the URL 

http : //www . nikhef . nl/~f orm/FORMdistr ibut ion/ packages/summer. 

The harmonic sums are related to the psi (digamma) function, which is the log- 
arithmic derivative of the Gamma function, and its derivatives, the polygamma 
functions ^ n \x), 

^(x)=^(x) = ^-\nT(x) = ^, (24) 
^ {n \x) = ^(x), (25) 



' Note that this is the convention used in [25]. Other authors denote a factor (— l) m 
by overlining the index. In this case a negative sign of the index indicates — conse- 
quently — a positive power of the corresponding variable. 
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via the following relations [29] : 



^<°>(n) = -7 + E \ = -7 + - 1), (26) 
fc=i ft 

^(™)(n) = (-l) m m! [-C(m + 1) + S m+1 (n - 1)] . (27) 

Using these relations the Gamma function can be expanded in terms of har- 
monic sums (cf. Ref. [30]) 

E|^±| = r(„) exp { £ £ - S,„(n - 1)} . (28) 



The expansion of the exponential function leads to the expansion formulae 
used in the following: 



r(i + e)r(n) 



r(l + e)r(n) 
T(n + e) 



= 1 +£51(71-1) (29) 
+ \e 2 [Sf(n-l)-S 2 (n-l)] 

+ | e 3 [Sf(n - 1) - SS^n - l)S 2 (n - 1) + 2S 3 (n - 1)] 
+0(^ 4 ), 

l-eSi(n-l) (30) 
+ ±5 2 [S 2 (n - 1) + S 2 (n - 1)] 

-J £ 3 [5 3 (n - 1) + 35'i(n - l)S 2 (n - 1) + 2S 3 (n - 1)] 
+0{e i ). 



3.2 Simplification of the sums 

The structure of the sum given in Eq. (15), 

oo oo l+m 

^ = EE E T(l,m,n), gem, (31) 

1=0 m=0 n =\l— m\ 
l+m+n=2g 

particularly the boundaries of the innermost sum, is very difficult to handle 
with symbolic manipulation programs. Therefore the first step is the simplifi- 
cation of this sum. 

For a given n = v we have to sum over pairs of indices (I, m), which satisfy 

\l -m\ < v < l + m. (32) 
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Neglecting the condition I + m + n = 2g at first, this is equivalent to 

m> I — v A m < I + v A m> v — I, (33) 
which is described by the rectangular area in the (/, m)-plane shown in Fig. 2. 

m = I + v 




Fig. 2. Area of summation for fixed n = v. 

Taking into account the condition / + m + n = 2g this means, that the sum 
ranges over all combinations of indices in this marked area whose sum is even. 

Let us in a next step consider pairs (7, m) with constant sum a — I + m, i.e. 
points on straight lines parallel to the second bisector. The sum starts with 
the line m = v — I [a = v) and reaches over the shaded rectangle in Fig. 2 
to infinity. Taking into account the symmetry of the terms it is sufficient to 
restrict oneself to the area below or over the first bisector. The contributing 
points on such a line can be described by the condition k < 21 < k + 2u, where 
u + k = a (cf. Fig. 2). 

If we parametrize k = 2k and k — 2k + 1, the addends read F(l, v + 2k — I, v) 
and F(l, v + 2k + 1 — I, v), respectively, and it is obvious, that in the second 
case the sum of the arguments can never be even, while the terms of the first 
kind automatically fulfill the condition. 

Therefore the original sum is transformed to a much simpler form: 

OO OO 71+K OO oo n 

F (l, n + 2K-l,n) = J2Y,J2 F (l + ^n + K-l,n). (34) 

n=l k=0 1=k n=l k=0 1=0 

If the expansions (29,30) are applied to the sum in Eq. (11), it is represented 
as a triple sum over rational expressions and products of such and harmonic 
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sums. In the following some algorithms are discussed to simplify such objects. 

3.3 Rational Expressions 

To sum expressions of the form 

5= f m 

^ (x + ai) (x + a 2 ) . . . (x + b x ) 2 (x + b 2 ) 2 ... (x + qi ) n (x + q 2 ) n 

(35) 

where P(x) is a polynomial in x, one first performs partional fractioning. This 
leads to simpler summands in which x appears only once in the denominator. 
The resulting sums can be evaluated easily in terms of zeta or psi functions 
by means of the formula [31] 

p 1 1/6 b \ 

^ (ax + 6)" a n V a a / 

= S^(^ +Q )-^ + " +1 ')- w 

If the upper limit is infinity, the second psi function is zero for n > 1. 
5.^ Algorithms for Nested Sums 

In Ref. [26] different algorithms have been published to evaluate nested sums 
using the properties of higher harmonic sums and more general structures 
like the so-called S- and Z-sums. For our purpose some special cases of the 
algorithms A and B of this work are needed. 

If an addend is not of the simple form S(a)/a m , but contains an offset in 
the denominator, which must not be symbolic, the summand can be reduced 
recursively 



Ski^) _y S kl (i) " 1 Sijz) | S kl (n) m 

~i(i + c) m ~i(i + c-l) m ~i(i + c-l) m i k {n + c) m ' 

This is a special case of Algorithm A. Of course, this is also valid for sums of 
depth greater than 2. 
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Algorithm B describes the recursive simplification of sums where products of 
harmonic sums appear. 



An expression of the form 



a=l 6=1 



can be rewritten by inserting the definition of Si. Then shifting variables, re- 
ordering of summations and partial fractioning yields to expressions which can 
be evaluated immediately and such of the same form as the original expression 
to which the procedure is applied again. Thus the sum is reduced recursively 
to simpler expressions and finally evaluated. 



3. 5 Products of Harmonic Sums 



In Ref. [32] algebraic relations between harmonic sums have been studied. For 
harmonic sums with identical arguments the permutation relation holds: 



Srn,n ^n,m S m S n -\- S m/ \ n , (39) 

where 

m A n := sign(m) sign(n) (\m\ + \n\). (40) 

This allows us to express the product of two harmonic sums in the following 

way: 

S m { x )S n (x) = S mjTl (x) + S njfn (x) — S m/ \ n (x) (41) 

For products with higher harmonic sums one has relations like 

Sa(x) ■ S bc (x) = S abc (x) + S bac (x) + S bca (x) - S aAb>c (x) - S b>aAc (x). (42) 
More relations of this kind can be found in [33] . 

3.6 Harmonic Sums with Double Argument 



In the following also sums with double argument Sk{2m) will appear. These 
cannot be handled with the algorithms described up to now. In Ref. [34] 
transformations of such objects to simpler sums have been discussed. From 
this work we need the relation 

oo i oo -i oo / I \i 

S ji~xr Sk{2j) = 2 " t_1 g -» Sk{l + 2x) + 2m_1 g -^ Ski * + 2x) (43) 
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which relates the sum over S(2k) to two sums over S(k), one of which is 
alternating. The opposite operation, the transformation of a sum with simple 
argument to one with double argument, can be done by means of a special 
case of the refinement formula [34] 

S m (n) = 2— 1 (S m (2n) + S_ m (2n)) . (44) 



Generally, we can write a sum of the form 

oo 

£/(2i) 

in the following way: 



(45) 



i=i 



£/(*)= E /(0 = 2 

igW ieven 



E/(0-E(-i)V(0 

Uew iew 



(46) 



The sum over all even numbers is replaced by the sum over all integers minus 
the sum over the odd numbers. This is achieved by subtracting an alternating 
sum and the division by 2. 



3.7 Multiple Infinite Sums 



Finally, we need some formulae to convert nested sums to multiple infinite 
sums or vice versa. Two important relations are given by 



oo I 



E E <f>( l > m ) = EE^ + m > m )' 

1=0 m=0 m=0 1=0 



(47) 



oo 21 



=1 m=l 



^^^ + m + l,2m + l)+^^^ + ffl + l,2m + 2). (48) 

m=0 1=0 rn=0 1=0 



Eq. (48) can be obtained by splitting the sum on the left-hand side into two 
sums, which is visualized in Fig. 3, and paramaterizing the summation vari- 
able m — 2/i and m = 2/x + 1, respectively. This leads to sums of the type 
J2m=o Sr^=m- These can be converted to a double infinite series by shifting the 
inner summation variable. 
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m 





im = 21 


• • • 




• • • / 




• • 













,m = 21 



m = 2l 




I 

Fig. 3. Splitting of a sum of type (48) in two sums 

If the argument of the appearing harmonic sums is a sum of two or more 
variables, the following relations are helpful: 



£ £ <f>(l, m) S(l + m) = J2Yl M n ~ 

/=0 rrt=0 n=0 Z=0 



(49) 



£5>(Z,m) S(l + 2m) 

1=0 m=0 

oo n oo n 

= Y,J2 - 2m - m ) S ( 2n ) + ^( 2n - 2m + 1, m) S(2n + 1). 

n=0 m=0 n=0 m=0 

(50) 



4 Convergence Acceleration 



If the sums cannot be completely simplified analytically, one can at least get 
a numerical result by summing up the terms with the computer. However, in 
our case the sum converges very slowly. In this Section we briefly sketch the 
method of nonlinear sequence transformations, by means of which one can 
improve the convergence of a series. For further information we refer to Ref. 
[35] and references cited therein. 

We consider a sequence {a n } and the corresponding infinite series 



S = J2a k 

k=0 



(51) 



with the partial sums 

S n = J2 ( 52 ) 

fc=0 

which themselves form a sequence {S n } that converges to the value S. The 
question is now, whether a different sequence {S' n } (or {a' n }) exists that con- 
verges to the same value S but has a better convergence behaviour. 



14 



The sequence {S' n } is said to converge faster than {S n } to the common limit 
S, if 

lim = o. (53) 

n ^°° s n — s 

In practice this means that S' N is a better approximation to the exact value 
of the limit than Sn, if the partial sum of order N is known. 

Convergence acceleration methods can be found with the help of model se- 
quences. This means, one investigates a particular class of sequences and con- 
structs a transformation, i.e. a prescription how to calculate the elements of a 
new sequence {S' n } out of the elements of {S n }, which is then exactly valid for 
the considered sequences. The prescriptions one has found in this way are often 
able to improve the convergence of different sequences of similar structure. 

One of the oldest and best known methods is Aitken's A 2 -process. The partial 
sums are transformed according to the rule 

c' c (AS n ) 2 (S n+ i — S n ) 2 

S n = s n - A2C = z n - — — — , n e JN . (54) 



A is the forward difference operator, which is defined as [35,31,27] 

Af(n) = f(n+l)-f(n). (55) 
Powers of A are to be understood as multiple application of A, e.g. 

A 2 f{n) = A [Af(n)\ = A [f(n + 1) - f{n)} 

= /(n + 2)-2/(n + l) + /(n). (56) 

Another method which can be found in the literature on numerical methods 
is Wynn's epsilon algorithm, which is defined by the recursion scheme 



6^ = 

et ] = S n (57) 



We follow the notation of Ref. [35]. The lower index indicates the transfor- 
mation order, the upper index labels the sequence elements. The difference 
operator acts only on the upper index. 
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Similar to the epsilon algorithm is Wynn's rho algorithm: 



P^ = S n (58) 

JP) _ (n+1) Xn+k+l - X n 

Pk+l-Pk-1 + _ (n7l) (n)' fc ' nG1N »' 

Pk ~Pk 



where the interpolation points Xk have to fulfill < x; t < Xj for i < j and the 
sequence {x n } diverges to infinity. A possible choice is x n — n + 1. 

An example, where the rho algorithm is very powerful, is the series defining 
C(2) = YhIi jz- When one sums the series up to the order % = 2000 one gets 
only 4 correct decimal digits. Application of the rho algorithm delivers more 
than 250 correct digits. 

For the Levin transformations remainder estimations have to be taken into 
account. The Levin transformation ist exact for sequences of the form 

k— i 

S n = S + u n J2 , C i a , v k,neli , (59) 

p (n + py 

where S is the limit of the sequence and (3 > is a parameter which can, 
in principle, be chosen arbitrarily and is usually set equal to 1. u n are the 
remainder estimates. They are, of course, generally not known exactly, because 
the actual limit is unknown. Therefore one uses assumptions for the u n which 
are customized to the given problem and which define different kinds of Levin 
t r ansf or mat ions . 

As an example, the ansatz 



LO n = , n e 1N , (60) 



*n+l 



yields Levin's v transformation 



MS ' \o) (/3+n+fc)fe-i a n+j a n+j+1 

v ( :\p, s n ) = ^ . (6i) 

V ( — 1 V ( k \ (/ 3 + n +J') fc 1 an+j-a n +j+l 
/=0 ^' (' 3 + n + k ) k ' 1 o, n+j a n+j+1 
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An enhancement of the epsilon algorithm is Brezinski 's theta algorithm: 




+ 



1 



k,n G 1N . 



(62) 



It is now suggestive to apply such a procedure, that improves the convergence 
of a series, succesively. In Ref. [35] algorithms are discussed which allow the 
calculation of high transformation orders efficiently and with low memory 
costs. 

The principle is to calculate for each new element of the series all possible ele- 
ments of the matrix d\. So one always gets the highest possible transformation 
order. Details can be found in Ref. [35], where also many more acceleration 
methods are discussed. 

For this work the algorithms mentioned above were programmed in FORTRAN. 
To achieve the desired numerical precision we used Bailey's library MPFUN 
[36,37], which allows calculations with arbitrary precision. 



5 Numerical Result for the Nonplanar Topology 

In this Section we sketch the application of the described methods to the 
nonplanar diagram iV and show how a numerical result is obtained. 

The sum we have to calculate reads (15) 



N s = {-As - Ae 2 + 32e 3 + 0(e 4 ))N s , 



(63) 



oo oo l+m 



N s = T,H E T(e,l,m,n,g), 



(64) 



1=0 m=0 n =\l-m\ 
l+m+n=2g 
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with the addend terms 



ne,l,m,n,g) = n2 _ i l e) ^ {1 _ e) 

x 1 T(g + 2-2e) 

(I + 1 - e) (m + 1 - e) (n + 1 - e) T(g + 2 - e) 
T(g-l+l-e)T(g-m + l-e)T(g-n + l-e) 

r(g-l + l)r(g-m + l)T{g-n + l) R V> l > m > n )> ^ 

where 
R(e, /, m, n) = 

4 

(2 + l + m + n-Qe) (n - 4 e) (1 + l + n-Ae) (2 + / + m + n - As) 

4 

+ (1 + T) (A + l + m + n) (2 + l + m + n-Qe) (1 + m - 2 7) + " ' ^ 

consists of rational expressions resulting from the radial integration of Eq. 
(13). 



5.1 Coefficient of e 1 



The coefficient of e 1 is a triple sum of rational expressions, where no harmonic 
sums occur. The sum can be evaluated analytically, which we sketch briefly in 
this subsection. The sum over / can be performed easily and we are left with a 
double sum. In the summands we get harmonic sums with arguments m and 
m + n. 

The resulting rational terms and terms with S(m) can now be summed over 
n. The terms containing S{m + n) can be cast into an expression that can 
be summed over m by means of a transformation of variables. After the sum- 
mation we are left with a simple sum with increased number of terms. The 
summands now consist of rational expressions, harmonic sums with simple 
argument Sk(i + n), such with double argument Sk(i + 2n) and products of 
the latter. 

Again, the rational expressions can be summed up analytically by means of 
Mathematica or SUMMER. Expressions containing harmonic sums with simple 
arguments can be summed with SUMMER. 

Expressions with harmonic sums with double argument can be transformed 
to terms calculable with SUMMER by means of formula (43). 
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A bit more complicated are terms where 2n appears in the denominator, like 

S 2 (n) 



S(2n + 1)«- 

By means of Eq. (44) the harmonic sum is converted to sums with double 
argument 

OO r\ 

E (2n+1)2 [S-2(2n) + S 2 (2n)}. (67) 
Then the sum is rewritten according to Eq. (46) 

tr^^[S-2(n) + S 2 {n)] (68) 

and can now be calculated. 
Finally, there are the products 

S k (n + i)Si(2m + j) (69) 
which can be evaluated with the help of Eq. (44) . 

The problem for the e 1 coefficient can thus be solved completely and the result 
reads 

iV 5 | £l = 68C 3 2 - 80C5 + 50C6, (70) 

which agrees with [10]. 
5.2 Coefficient of e 2 



In the expression for the e 2 coefficient the terms in the triple sum contain the 
harmonic sums Si(n — I), Si(m), Si(l) and Si(m + n). All expressions can be 
summed over a variable which does not appear as an argument in S. 

The result is then a double sum 

££F(m,n), (71) 

m=l n=l 

where terms with harmonic sums of the following types occur: 

Si (m) , Si(m) Sj (m) , S^m) Sj (n) , 

Si(m + n), Si(m + n) Sj(m + n), Si(m) Sj(m + n). 
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In the case of Si(m) and Si(m) Sj(m) we sum now over the variable n which is 
not the argument of S. Expressions containing Si(m + n), Si(m + n) Sj(m + n) 
can be handled by means of variable transformations as described in Section 
3.7. The terms with Si(m) Sj(n), S'j(m) Sj(m + n) can in the most cases be 
cast into a form to which algorithm B of Ref. [26] can be applied. 

Thus the second summation yields a simple sum. Now the addends contain 
simple harmonic sums and products of two or three harmonic sums, some of 
which have double arguments: 

Si(m), Si(m) Sj(m), S i (m)S j (2m), S^m) Sj(m) S k (2m). 

These terms can in principle be calculated by means of the described algo- 
rithms. 

However, in some of the terms there are certain combinations of denomina- 
tors and harmonic sums in which the summation variable is doubled. These 
can in some cases be be handled with the known algorithms, however, not 
in the general case. Therefore the problem of the order e 2 cannot be solved 
analytically at present. To achieve this, more examinations in the subject of 
the summations are necessary. 



5.3 Numerical evaluation 

Since up to now an analytical solution is not possible with our method we 
try to find a numerical solution. To do this, we simplify the terms as far as 
possible and calculate the sums numerically. 

For the numerical calculation the terms have been simplified with Mathematica 
and output in FORTRAN syntax. This code has been optimized by means of 
Maple and built in a FORTRAN program. To speed up the calculation this pro- 
gram was parallelized by means of MP I [38,39]. Finally, the program was trans- 
lated with Bailey's program TRANSMP [37] to allow arbitrary numerical preci- 
sion by means of the MPFUN library [36]. The programs store all partial sums. 
To the sequence of partial sums we apply different convergence acceleration 
methods. 

A different approach to a numerical solution, is the following: We begin with 
the representation of the sum (34) 

oo oo n 

J2F(l + m,n + m-l,n). (72) 

n=l m=0 1=0 
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For a fixed n — v the expression 

V oo 

J2 J2 F(l + m,v + m-l,is) (73) 

1=0 m=0 

describes v + 1 infinite sums 

oo 

S v ,i = J2 F(l + m,is + m-l,is). (74) 

m=0 



Sums of that kind can easily be calculated with SUMMER. So the obvious pro- 
cedure is to evaluate all contributions S n j for all values of n up to an upper 
limit n max . 

Thus we get a sequence of expressions consisting of rational numbers and zeta 
functions. These can be evaluated numerically to achieve a sequence of decimal 
numbers {S n }, where 

lim S n = N s \s 2 - 

To this we again apply the convergence acceleration algorithms. The first 
elements of this sequence are 

5 = 56 ( 2 - 56 Cs + 32 ( 2 Cs - 80 Cf - 100 ( 4 + 204 ( 3 C 4 

+64 C 5 + 144 C 2 Cs - 488 Ce + 308 C 7 , 
1220 848 644 

51 = — ^-C 2 - — C 3 + 16OC2C3 - sod - — C 4 + 204C 3 C 4 

3232 

+320 C 5 + 144 C 2 C 5 - 488 Ce + 308 C 7 + — , 

52 = - ^ ~ 2 ~W Cs + 1 6 8 C2 Cs ~ 8 Cl " IT C4 + 2 04 Cs C4 

64279 

+396 C 5 + 144 C 2 C 5 - 488 Ce + 308 (7 + 
and the 500th element is (with rounded coefficients) 



S500 = 109.191 - 65.7377 C 2 - 114.167 C 3 + 210.296 C2C3 ~ 80. C3 - 741.581 ( 4 
+204. Cs C4 + 500.591 Cs + 144. ( 2 Cs - 488. Ce + 308. ( 7 
^ 205.62576485415736523. 



This procedure can be implemented easily in FORM or its parallelized version 
PARF0RM [40,41]. 

We show in Table 1 the results of the FORTRAN-Program. The entry in the line 
labeled 'Sum' is the result we get by simply summing up the terms, where the 
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upper limit of the outermost summation was set to 15,000. The other lines 
contain the results obtained by application of different algorithms. 



Method 


Value for Ng \ £ 2 


Sum 


205.625 765 027 123 610 885 997 875 


Theta 


205.625 765 027 124 189 154 823 016 


Epsilon 


205.625 765 027 124 195 743 663 390 


Aitken 


205.625 765 027 124 195 744 291 025 


Rho 


205.625 765 027 124 195 744 492 088 


Levin v 


205.625 765 027 124 196 134 222 058 



Table 1 

Numerical results for the e 2 coefficient of N 

One can see, that the results of the Aitken and the Rho method agree in 21 
decimal digits. These algorithms delivered the best result * for the numerical 
calculations of the planar topology, where the exact result was known and 
could be used to check the results. Since the structure of the terms is similar, 
one can expect, that they are also suitable for the nonplanar diagram. The 
result of the Epsilon procedure agrees in 20 digits with the former. All results 
agree in 16 digits. In a very pessimistic interpretation at least these digits can 
be considered reliable. 

The results of the numerical evaluation by means of the FORM-method agree 
with these numbers in 14 or 20 digits, depending on the acceleration method 
used. 

5.4 Integer Relation Detection 

Finally, we want to get an idea of the analytical result by means of "exper- 
imental mathematics". Looking at the known result for the planar diagram, 
we assume, that the result for N be also a linear combination of zeta values 
with integer coefficients, which have to be determined. This problem is known 
as integer relation detection [42]. 

In general a vector x = (xi, x 2 , ■■■,x n ) with real or complex elements is said to 
possess an integer relation, if integers exist such that at least one of them 
is different from zero and 

a l x l + a 2 x 2 + . . . + a n x n = (75) 

* Summing up 2000 terms we get 3 correct digits for P. Application of the rho 
algorithm yields 13 correct digits in this case. 
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An algorithm which can solve this problem is the PSLQ algorithm [43,44,45]. 

The input for this algorithm are in our case the numerically calculated number 
and numerical values for the constants, i.e. the zeta values. Then the program 
checks, whether there is an integer relation between the numbers within the 
desired accuracy. We used an implementation of this method by O. Veretin 
[46]. 

Assuming that the transcendentals are the same as for the planar diagram, as 
it is the case in the lower orders in e, we find the following result for the 2nd 
order coefficient of the e expansion of N: 

N s \ £ 2 = -272 Cl + 204 ( 3 ( 4 + 80( 5 - 200( 6 + 450( 7 . (76) 

This agrees with the analytical result which was found in the meantime [19,20]. 

Comparison of this result with Table 1 shows, that all the digits printed for 
the rho algorithm are correct. This proves the power of this algorithm for our 
problem and justifies the confidence in the convergence acceleration techniques 
retrospectively. 



6 Analytical Result for the Planar Topology 

In this Section we show that the coefficient of e 2 for the planar diagram can 
be found analytically by means of the described method. The steps we take 
are the same as described in Section 5. The structure of the sums and of the 
terms are similar and so is the procedure of the calculation. The difference 
is, the terms that cause problems in the case of N, because we do not know 
algorithms for their evaluation, do not appear in the case of P. 

After application of the algorithms we have the sum in the form 

oo 

Ps\ £ 2 = P + Pi + Y, p 2(a), (77) 

o=l 

where P is the contribution for n = 0, 

P = 56 Ca - 164 C 4 + 56 Cs - 440 Ce + 320C 7 

+ 48 C 2 Cs + 180 C 3 C4 + 120 C 2 C 5 - 80 C 3 2 , (78) 
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and 



17686997 166771 86459 , 

P l = (2 + Ca - 192 Ci C3 + H52 C2 Ca 

1 17496 324 S2 162 S3 U S3 S2S3 

2702 

- I6C1C2C3 - 592 C| + 128 Ci C 3 2 - — C4 + 552 Ci C4 

394 6770 

- 1032 C3C4 - —C 5 - 656CiC 5 - 1240 C 2 C 5 - ^-Ce 



+ 770 CiCe + 2248 C 7 , 



528 Si (a) 480Si(a) 48Si(a) 2 48Si(a)S 2 (a) 6Si 3 (a) 
P2[a) ~ ~^ ? a* + ~JT^ 

8 S'i(a) 2 S 2 (a) 32 S^a) S 2 (a) S 3 (a) _ 24 S^a) S 1;3 (a) 



(3 + ay 1 + a (3 + a)' 



P 2 (a) consists only of terms that can be calculated with SUMMER. Doing this 
we arrive at the known result 

P s \ £ 2 = -176 C 3 2 + 132 C 3 C 4 + 80 C 5 - 200 Ce + 317 Cr- (79) 



7 Conclusion 



We described a method for the calculation of massless, dimensionally regu- 
larized Feynman integrals in their e-expansion. It is based on the expansion 
of the integrand in terms of Gegenbauer polynomials. This yields sums over 
Gamma functions which can be expanded in harmonic sums. As a consequence 
one arrives at sums over rational functions and products of such with harmonic 
sums which can be simplified by the introduced summation algorithms and 
then be calculated numerically or analytically. 

The procedure in principle is independent of the Gegenbauer polynomial tech- 
nique. It can be useful wherever there are sums of the described form. The 
method should be expandable to more complicated problems, e.g. higher or- 
ders of the diagrams under consideration or more complex diagrams. 

In the case of the nonplanar diagram we could not get to a complete analytical 
solution. To achieve this, further knowledge in the summation techniques is 
required. 

For this work we combined the abilities of different computer programs, in 
particular the pattern matching facilities of Mathematica, the code genera- 
tion functions of Maple, and FORM's power of managing very large expressions. 
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It is the idea to implement the method in one programming language to make 
it stand-alone. This could either be done in FORM or in C++. The summation al- 
gorithms of Ref. [26] have already been implemented in both, the FORM version 
has recently been published [50]. In [51] also a Mathematica implementation 
was introduced. 



8 Summary of the results 

The coefficients of the e expansion of the nonplanar three-loop-diagram N 
could be calculated numerically to 20 digits 

N s = 20.738555102867398527 + 66. 168906981239990785 e 

+ 205.62576502712419574 e 2 

The analytical result for the 0(e) term given in [10] could be reproduced, 

N s = 20 C 5 + (68 C| - 80 C 5 + 50 Ce) e. (80) 

Application of the PSLQ algorithm to the 0(e 2 ) term gives 

N s \ £ 2 = -272 C| + 204 C 3 C 4 + 80Cs - 200Ce + 450C 7 , (81) 
which agrees with [19,20]. 

Also the analytical result up to order e 2 for the planar diagram was verified. 

P s y = -176 C 3 2 + 132 Cs C4 + 80 Cs - 200 Ce + 317 ( 7 . (82) 
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A Properties of the Gegenbauer Polynomials 

The Gegenbauer polynomials, also called ultraspherical polynomials, are the 
coefficients in the power series exansion of the generating function 

1 oo 

They are orthogonal on the interval [—1, +1] with the weight function 

w(x) = (1 -x 2 Y x ~^, (A.2) 



Special cases are the Chebyshew polynomials of first or second kind with 
A = 0, 1, respectively, and the Legendre polynomials with \ —\. 

Information about these polynomials can be found in the books by Erdelyi 
et al. [47], Vilenkin [48] or Szego [49]. The notation in this article follows 
Chetyrkin et al. [12]. 

From the orthogonality relation (A. 3) follows a relation for unit vectors x, y, z, 
which is used to reduce angular integrals in Section 2: [12,47] 

J dy C n A (x ■ y) C^(y • z) = -A_ 5 n , m C n A (x • z) (A.4) 



Values for special parameters are: 

C^(x) = 1, (A.5) 



T(n + 2X) 
n{ ' n!T(2A) : 



(A.6) 
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Q A (o) = 





HTW m - n 
m!T(A) ' 1,1 ~ 2 



n odd 



n even 



(A.7) 



Another basic feature of the Gegenbauer polynomials used in this work is the 
expansion of a propagator: 



where 



1 



1 



(x! — x 2 ) 2A max(r 1 ,r 2 ) ; 



oo / \ n/2 

E^(X!-X 2 ) ( n 



n=0 



\r 2 



x 



n\ (ri r 2 

— ) = mm — , — 

r 2 l \r 2 n 



r 1 /r 2 ifr! < r 2 
r 2 /r\ ifr 2 < r\ 



(A.8) 
(A.9) 

(A.10) 



The exponential function can be expanded by means of 

oo 

e 2, P x = r(A) j2 F (n + A) C*(A ■ p) (p 2 r)«/ 2 j x+n (p 2 r), 



n=0 



where j a (z) is related to the Bessel function J a (z) by 
and the following equation holds: 



(A.H) 



(A.12) 



Rea>2Re6 + i 



(A.13) 



Finally, a formula is needed which allows to rewrite the product of two Gegen- 
bauer polynomials as sum of single Gegenbauer polynomials. 



l+m 



Cf{x)C x m {x)= £ D x (l,m;n)C*(x), 



(A.14) 



n=\l—m\ 
l+m+n=2g,g£Z 



D\(l,m;n) 



n\ (n + \)T(g + 2\) 
T 2 (\)T(g + \ + l)T(n + 2\) 

T(g-l + A) T(g - m + A) T(g - n + A) 
X T{g-l + \)T{g-m + \)T{g-n + \) 



. (A.15) 
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